#define _USE_MATH_DEFINES
#include <cmath>
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <string>
#include <vector>
#include <sstream>
#include <ctime>

#include "srandInitializer.hpp"
#include "randomNumberGenerators.hpp"
#include "matrix.hpp"
#include "matrixdecomposition.hpp"
#include "timeSequence.hpp"
#include "varianceStructure.hpp"
#include "SimpleMultiDimensionGaussian.hpp"
#include "LMMAdvancer.hpp"
#include "LMMPricer.hpp"
#include "newmat10/newmat.h"
#include "newmat10/include.h" 


using namespace TermStructure;
using namespace std;

template<typename T>
std::ostream &outputTo(std::ostream &stream, T &m)
{
	for(int i = 0; i < m.GetLength(0); i++)
	{
		for(int j = 0; j < m.GetLength(1); j++)
		{
			stream<<m(i,j)<<",";
		}
		stream<<endl;
	}
	return stream;
}
int main(int argc, char const *argv[])
{
	clock_t start, stop;
	//Variance and time Sequence
	TermStructure::TimeSequence<double> times2;
	TermStructure::VarianceStructure<double> times1;
	//Put times into time sequence
	for (int i=1;i<=6;i++)
	{
		times2.PutTime(0.25*i);
	}
	//Set time sequence to variance
	times1.SetTimeSequence(times2);
	//Set Variance Parameters (according to the doc)
	times1.d=0.05;
	times1.a=0.09;
	times1.b=0.44;
	times1.c=0.11;
	times1.beta2=0.1;

	//Create the Advancer
	LMMAdvancer<
		VarianceStructure<double>,
		SimpleMultiDimensionGaussian,
		double> la;
	la.InitiateForwardRateSequence();
	//Set variance and time sequence
	la.Set(times1, times2);
	//Set current time to zero;
	la.CurrentTime = 0.0;
	//Set the time step to 0.25
	la.TimeStep = 0.5;
	//Tell advancer to pre-compute covariance matrix and set it to Gaussian Generator, therefore do not need to re-decompose the matrix again
	la.PreAdvance();
	//Set current forward rates to 0.05
	for(int i = 0; i < 7; i++) 
	{
		la.CurrentSequence()[i] = 0.05;
	}
	//Now we compute the expected value of forward rate from 0.25~1.5 years
	double t = 0;
	int N = 1;
	
	for(int i = 0; i < N; i++)
	{
		la.SingleAdvance();
		la.NextSequence().EstablishDiscountFactors();
		t += la.NextSequence().SwapRate(1, 6)/N;
	}
	cout<<t<<endl;

	//Now we compute the expected value of the caplet(0,0.25,0.5)@0.05, with a theoretical value of 0.000408154
	double c = 0;
	for(int i = 0; i < N; i++)
	{
		la.SingleAdvance();
		if(la.NextSequence()[2] >= 0.05) {
			la.NextSequence().EstablishDiscountFactors();
			c += times2.TimeDifference(2)*(la.NextSequence()[2] - 0.06>0?la.NextSequence()[2] - 0.05:0)/(1+la.NextSequence()[2]*times2.TimeDifference(2))/N/(la.NextSequence().DiscountFactorAt(6,1));
		}
	}
	la.CurrentSequence().EstablishDiscountFactors();
	cout<<c * la.CurrentSequence().DiscountFactorAt(6)<<endl;
	//	return 0;

	//Now we compute the expected value of the caplet(0,0.25,0.5)@0.05, with a theoretical value of 0.000408154


	TermStructure::SimpleMatrix<double> a(5,5);
	TermStructure::SimpleMatrix<double> u(5,5);
	TermStructure::SimpleMatrix<double> d(5,5);
	TermStructure::SimpleMatrix<double> v(5,5);
	for(int i = 0; i < 5; i++)
	{
		for(int j = 0; j < 5; j++)
		{
			a(i,j) = 1.0/(i+j+1);
		}
	}
	SingularValueDecomposeTo(a,u,d,v);	

	for(int i = 0; i <5; i++)
	{
		for(int j = 0; j <5; j++)
		{
			std::cout<<u(i,j)<<" ";
		}
		cout<<endl;
	}
//	return 0;

	SimpleMatrix<double> var(5,5);
	for(int i = 0; i < 5; i++)
	{
		for(int j = 0; j < 5; j++)
		{
			if(i==j) var(i,j) = 1;
			else var(i,j) = 0.9;
		}
	}
	TermStructure::PCABasedMultiDimensionGaussian<SimpleMatrix<double>,double> pca;
	pca.SetVarianceMatrix(var);
	ofstream out("tester.csv");
	for(int i = 0; i < 10000; i++)
	{
		pca.NextGroup();
		for(int j = 0; j < 5; j++)
		{
			out<<pca(j,0)<<",";
		}
		out<<endl;
	}
	out.close();




	ifstream inFile ("test.csv");
    vector<double> tmp;
	string line;
    int linenum = 0;
    while (getline (inFile, line))
    {
		
        linenum++;
       // cout << "\nLine #" << linenum << ":" << endl;
        istringstream linestream(line);
        string item;
        int itemnum = 0;
        while (getline (linestream, item, ','))
        {
            itemnum++;
        //    cout << "Item #" << itemnum << ": " << item << endl;
			if(itemnum == 2) {tmp.push_back(atof(item.c_str())); };
		}
    }
	

	TimeSequence<double> ts;
	for(int i = 1; i <tmp.size() ; i++)
	{
		ts.PutTime(i*0.25);
	}

	ForwardRateSequence<double> frs;
	frs.SetTimeSequence(ts);
	for(int i = 0; i < ts.GetSize(); i++)
	{
		frs[i] = tmp[i];
	}

	VarianceStructure<double> vs;

	vs.d=0.3135;
	vs.a=0;
	vs.b=0.8146;
	vs.c=0.4814;
	vs.beta1=0.4465;
	vs.beta2=0.8791;

	
	LMMPricer<VarianceStructure<double>,SimpleMultiDimensionGaussian,double> pricer; //PCABasedMultiDimensionGaussian
	pricer.SetLMMAdvancer(ts, vs);
	ofstream out1("result.csv");
	pricer.InitialForwardRateSequence() = frs;
	Swap<double> swap1(ts,20, 10);
	//cout<< pricer.PriceEuropeanSwaption(swap1, 0, 10,10000);
/*	
	for (int ss =1;ss<160;ss++)
	{
		Swap<double> swap(ts,160, ss);
		//Swap<double> swap(ts, pricer.InitialForwardRateSequence().GetSize(), ss);
	
		out1 << pricer.PriceEuropeanSwaption(swap, 0, ss,10000)<<",";
		cout << ss << " : " << endl;
	//cout<<"European Swaption:"<<pricer.PriceEuropeanSwaption(swap, 0, ss)<<endl;
	}

*/
	Swap<double> swap(ts,20, 4);
	TimeSequence<double> chooser;
	for (int i =1;i<4;i++)
	{
		chooser.PutTime(i*0.25);
	}
	
 	 start = clock(); 
	// out1 << "European Swaption : " << pricer.PriceEuropeanSwaption(swap, 0, 4,5000) << endl;
	 stop = clock(); 
	 out1 << double(stop - start) / CLOCKS_PER_SEC <<endl;
	 cout << "done" <<endl ;  


	 start = clock(); 
	// out1 << "Trigger Swaption : " << pricer.PriceTriggerSwap(swap,0,chooser,5000)<<endl;
	 stop = clock(); 
	 out1 << double(stop - start) / CLOCKS_PER_SEC <<endl;
	 cout << "done" <<endl ;  
	
	 start = clock(); 
//	 out1 << "Ratchetcaps :" << pricer.PriceRatchetcaps(5*4,5000) << endl;
	 stop = clock(); 
	 out1 << double(stop - start) / CLOCKS_PER_SEC <<endl;
	 cout << "done" <<endl ;  


	 start = clock(); 
//	 out1 << "Bermudan Swaption : " <<pricer.PriceBermudanSwaption(swap)<<endl;
	 stop = clock(); 
	 out1 << double(stop - start) / CLOCKS_PER_SEC <<endl;
	 cout << "done" <<endl ;  
	 

	 start = clock(); 
//	 cout << "Const Maturity Swap : " << pricer.PriceConstMaturitySwap(7*4,swap,2*4,0,1000)<<endl;
	 stop = clock(); 
	 out1 << double(stop - start) / CLOCKS_PER_SEC <<endl;
	 cout << "done" <<endl ;  

	 start = clock(); 
//	 cout << "Price Caption : " << pricer.PriceCaption(10*4,0.005,0.000001,5*4,5000)<<endl;
	 stop = clock(); 
	 out1 << double(stop - start) / CLOCKS_PER_SEC <<endl;
	 cout << "done" <<endl ;  

	out1.close();
 

	 start = clock(); 
	 double cms_30 = pricer.PriceConstMaturitySwap(7*4,swap,2*4,0,5000);
	 double cms_10 = pricer.PriceConstMaturitySwap(7*4,swap,30*4,0,5000);
	stop = clock(); 
	 cout << "Const Maturity Swap : " << cms_30 - cms_10<<endl;
	 
	 cout << double(stop - start) / CLOCKS_PER_SEC <<endl;

	 cout << "done" <<endl ;  

	/*


	TimeSequence<double> chooser;
	for (int i =1;i<4;i++)
	{
		chooser.PutTime(i*0.25);
	}
*/
//	Swap<double> swap(ts);

//	cout<<"Trigger Swaption:"<<pricer.PriceTriggerSwap(swap, 0.05, chooser,5000)<<endl;
//	cout<<"Const Maturity Swap"<<pricer.PriceConstMaturitySwap(10,swap1, 22,2,0,1000)<<endl;
//	cout<<"Ratchet caps"<<pricer.PriceRatchetcaps(10,1000)<<endl;
//	cout<<"caplet:"<<pricer.PriceCaplet(0,0.25,0.5,0.05,100000)<<endl;
//	cout<<"swap:"<<pricer.PriceSwap(swap,0,40,4,0.05,100000)<<endl; //swap(0;1,2.5;K=0)

	return 0;

}